---
title: "CARRS_ByEducation"
author: ""
date: ""
output: html_document
---


``` {r BasicStuff}
Sys.setenv(LANGUAGE = "en")
options(scipen = 999)
```

```{r , setup, include=FALSE, eval = TRUE}
knitr::opts_chunk$set(
   comment = '', warning=FALSE,  message=FALSE, results= FALSE , fig.width = 9,fig.height = 7, echo=FALSE, root.dir =""
)

 knitr::opts_knit$set(
   root.dir = ""    
 )
```


``` {r Libraries, warning = FALSE, message = FALSE, include = FALSE}
library(dplyr)
library(tidyr)
library(data.table)
library(haven)
library(tidyverse)
library(ggplot2)
library(rdlocrand)
library(rddensity)
library(rdrobust)
library(openxlsx)
library(ggtext)
library(extrafont)
library(viridis)
library(scales)
library(Synth)
library(broom)
library(fixest)
library(MASS)
library(scales)
library(extrafont)
library(modelsummary)
library(kableExtra)
library(ggpubr)
loadfonts()

```


# Loading dataset

``` {r RDdata, eval = TRUE}
# data for outcome: drug treatment
  CARRSfull<-as.data.table(read.csv("carrs_full.csv"))
  CARRSfull <- na.omit(CARRSfull, cols=(c("female","hypdrug","edu")))  
 
# data for outcome: diagnosis
  CARRSsel<-as.data.table(read.csv("carrs_selected.csv"))
  CARRSsel <- na.omit(CARRSsel, cols=(c("female","hypdiag","edu")))
  
# data for outcomes: systolic and diastolic blood pressure
  CARRSbp <-as.data.table(read.csv("carrs_bp.csv"))
  CARRSbp <- na.omit(CARRSbp, cols=(c("sysdiff","diasdiff","edu")))
```

# RDD function - binary outcomes

``` {r RDfunction, eval = TRUE }
RDD_educ <- function(Y) {
  out =  rdrobust(y = Y, X, c = 0, cluster = finaldata$pid, covs = cbind(finaldata$w01, finaldata$w23, finaldata$w45),
                  kernel = "triangular", p = 1, rho = 1, all = T)
  ret <- data.frame(
    res = c(formatC(round(100*out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(100*out$ci[3,1],2),2,format="f"), " - ", formatC(round(100*out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2]))
  )
  rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")

  # needed for the figures
  estimate <- out$coef[1,1]
  CIl <- out$ci[3,1]
  CIu <- out$ci[3,2]
  
  return(list(ret=ret, estimate=estimate, CIl=CIl, CIu=CIu))
}
```

# Outcome: diagnosed in the past 12 months

``` {r DiagEdu, eval = TRUE }
# Those less than secondary school
  finaldata <- CARRSsel[edu == 0,]

  X = finaldata$runningvar_edu0
  outcome = finaldata$hypdiag
  outLow <- RDD_educ(Y = outcome)
  
  resL <- outLow$ret
  
  resL <- resL %>%
    rename(Low = res)%>%
    add_rownames

  figL <- c("Low", "1 total",(outLow$estimate)*100, (outLow$CIl)*100, (outLow$CIu)*100)

# Those with high school
  finaldata <- CARRSsel[edu == 1,]

  X = finaldata$runningvar_edu1
  outcome = finaldata$hypdiag
  
  outMiddle <- RDD_educ(Y = outcome)

  resM <- outMiddle$ret
  
  resM <- resM %>%
    rename(Middle = res)

  figM <- c("Middle", "1 total",(outMiddle$estimate)*100, (outMiddle$CIl)*100, (outMiddle$CIu)*100)

# Those with more than high school 
  finaldata <- CARRSsel[edu == 2,]

  X = finaldata$runningvar_edu2
  outcome = finaldata$hypdiag
  
  outHigh <- RDD_educ(Y = outcome)
  
  resH <- outHigh$ret
  
  resH <- resH %>%
    rename(High = res)

  figH <- c("High", "1 total",(outHigh$estimate)*100, (outHigh$CIl)*100, (outHigh$CIu)*100)

  EduDiag <- cbind(resL,resM,resH)
  figEdu <- rbind(figL, figM, figH)
```

``` {r DiagEduSex, eval = TRUE }
# Those with primary school or less
  finaldata <- CARRSsel[edu == 0 & female == 0,]

  X = finaldata$runningvar_male_edu0
  outcome = finaldata$hypdiag
  
  outLowM <- RDD_educ(Y = outcome)
  
  resLM <- outLowM$ret
  
  resLM <- resLM %>%
    rename(Low = res)%>%
    add_rownames
  
  figLM <- c("Low", "2 men",(outLowM$estimate)*100, (outLowM$CIl)*100, (outLowM$CIu)*100)  

  
  finaldata <- CARRSsel[edu == 0 & female == 1,]

  X = finaldata$runningvar_female_edu0
  outcome = finaldata$hypdiag
  
  outLowF <- RDD_educ(Y = outcome)
  
  resLF <- outLowF$ret
  
  resLF <- resLF %>%
    rename(Low = res)%>%
    add_rownames

  figLF <- c("Low", "3 women",(outLowF$estimate)*100, (outLowF$CIl)*100, (outLowF$CIu)*100)  

# Those with secondary school
  finaldata <- CARRSsel[edu == 1 & female == 0,]

  X = finaldata$runningvar_male_edu1
  outcome = finaldata$hypdiag
  
  outMiddleM <- RDD_educ(Y = outcome)

  resMM <- outMiddleM$ret
  
  resMM <- resMM %>%
    rename(Middle = res)
 
  figMM <- c("Middle", "2 men",(outMiddleM$estimate)*100, (outMiddleM$CIl)*100, (outMiddleM$CIu)*100)
  
   
  finaldata <- CARRSsel[edu == 1 & female == 1,]

  X = finaldata$runningvar_female_edu1
  outcome = finaldata$hypdiag
  
  outMiddleF <- RDD_educ(Y = outcome)

  resMF <- outMiddleF$ret
  
  resMF <- resMF %>%
    rename(Middle = res)

  figMF <- c("Middle", "3 women",(outMiddleF$estimate)*100, (outMiddleF$CIl)*100, (outMiddleF$CIu)*100)

  # Those with more than high school 
  finaldata <- CARRSsel[edu == 2 & female == 0,]

  X = finaldata$runningvar_male_edu2
  outcome = finaldata$hypdiag
  
  outHighM <- RDD_educ(Y = outcome)
  
  resHM <- outHighM$ret
  
  resHM <- resHM %>%
    rename(High = res)

  figHM <- c("High", "2 men",(outHighM$estimate)*100, (outHighM$CIl)*100, (outHighM$CIu)*100)

  
  finaldata <- CARRSsel[edu == 2 & female == 1,]

  X = finaldata$runningvar_female_edu2
  outcome = finaldata$hypdiag
  
  outHighF <- RDD_educ(Y = outcome)
  
  resHF <- outHighF$ret
  
  resHF <- resHF %>%
    rename(High = res)
  
  figHF <- c("High", "3 women",(outHighF$estimate)*100, (outHighF$CIl)*100, (outHighF$CIu)*100)
  
  EduMDiag <- cbind(resLM, resMM, resHM)
  EduFDiag <- cbind(resLF, resMF, resHF)
  
  figSexEdu <- rbind(figLM, figMM, figHM, figLF, figMF, figHF)
  

  # dataset for figures with all results
  DiagnosisResultsEdu <- rbind(figEdu,figSexEdu)
  DiagnosisResultsEdu <- as.data.frame(DiagnosisResultsEdu)
  DiagnosisResultsEdu <- rename(DiagnosisResultsEdu, edu = V1, bp.sample = V2, estimate = V3, CIl = V4, CIu = V5)
  DiagnosisResultsEdu <- transform(DiagnosisResultsEdu, estimate = as.numeric(estimate), CIl = as.numeric(CIl), CIu = as.numeric(CIu))
  DiagnosisResultsEdu$bp.sample = factor(DiagnosisResultsEdu$bp.sample, levels = c("3 women", "2 men", "1 total"))
  DiagnosisResultsEdu$edu = factor(DiagnosisResultsEdu$edu, levels = c("High", "Middle", "Low"))
```

``` {r DiagAge, eval = TRUE }
# 30-39 years
  finaldata <- CARRSsel[agecat == 0,]

  X = finaldata$runningvar_age0
  outcome = finaldata$hypdiag
  outYoung <- RDD_educ(Y = outcome)
  
  resY <- outYoung$ret
  
  resY <- resY %>%
    rename(Young = res)%>%
    add_rownames 

  figY <- c("30-39 years", "1 total",(outYoung$estimate)*100, (outYoung$CIl)*100, (outYoung$CIu)*100)

# 40+ years
  finaldata <- CARRSsel[agecat == 1,]

  X = finaldata$runningvar_age1
  outcome = finaldata$hypdiag
  
  outOld <- RDD_educ(Y = outcome)

  resO <- outOld$ret
  
  resO <- resO %>%
    rename(Old = res)
  
  figO <- c("40+ years", "1 total",(outOld$estimate)*100, (outOld$CIl)*100, (outOld$CIu)*100)

  AgeDiag <- cbind(resY, resO)
  figAge <- rbind(figY, figO)
```

``` {r DiagAgeSex, eval = TRUE }
# 30-39 years
finaldata <- CARRSsel[agecat == 0 & female == 0,]

  X = finaldata$runningvar_male_age0
  outcome = finaldata$hypdiag
  
  outYoungM <- RDD_educ(Y = outcome)
  
  resYM <- outYoungM$ret
  
  resYM <- resYM %>%
    rename(Young = res)%>%
    add_rownames
  
  figYM <- c("30-39 years", "2 men",(outYoungM$estimate)*100, (outYoungM$CIl)*100, (outYoungM$CIu)*100)  

finaldata <- CARRSsel[agecat == 0 & female == 1,]

  X = finaldata$runningvar_female_age0
  outcome = finaldata$hypdiag
  
  outYoungF <- RDD_educ(Y = outcome)
  
  resYF <- outYoungF$ret
  
  resYF <- resYF %>%
    rename(Young = res)%>%
    add_rownames 
  
  figYF <- c("30-39 years", "3 women",(outYoungF$estimate)*100, (outYoungF$CIl)*100, (outYoungF$CIu)*100)  

# 40+ years
  finaldata <- CARRSsel[agecat == 1 & female == 0,]

  X = finaldata$runningvar_male_age1
  outcome = finaldata$hypdiag
  
  outOldM <- RDD_educ(Y = outcome)

  resOM <- outOldM$ret
  
  resOM <- resOM %>%
    rename(Old = res)

  figOM <- c("40+ years", "2 men",(outOldM$estimate)*100, (outOldM$CIl)*100, (outOldM$CIu)*100)
  
  finaldata <- CARRSsel[agecat == 1 & female == 1,]

  X = finaldata$runningvar_female_age1
  outcome = finaldata$hypdiag
  
  outOldF <- RDD_educ(Y = outcome)

  resOF <- outOldF$ret
  
  resOF <- resOF %>%
    rename(Old = res)

  figOF <- c("40+ years", "3 women",(outOldF$estimate)*100, (outOldF$CIl)*100, (outOldF$CIu)*100)
  
  AgeMDiag <- cbind(resYM, resOM)
  AgeFDiag <- cbind(resYF, resOF)

  figSexAge <- rbind(figYM, figOM, figYF, figOF)

  # dataset for figures with all results
  DiagnosisResultsAge <- rbind(figAge,figSexAge)
  DiagnosisResultsAge <- as.data.frame(DiagnosisResultsAge)
  DiagnosisResultsAge <- rename(DiagnosisResultsAge, age = V1, bp.sample = V2, estimate = V3, CIl = V4, CIu = V5)
  DiagnosisResultsAge <- transform(DiagnosisResultsAge, estimate = as.numeric(estimate), CIl = as.numeric(CIl), CIu = as.numeric(CIu))
  DiagnosisResultsAge$bp.sample = factor(DiagnosisResultsAge$bp.sample, levels = c("3 women", "2 men", "1 total"))
  DiagnosisResultsAge$age = factor(DiagnosisResultsAge$age, levels = c("40+ years","30-39 years"))
```

# Outcome: currently in treatment

``` {r TreatEdu, eval = TRUE }
# Those less than secondary school
  finaldata <- CARRSfull[edu == 0,]

  X = finaldata$runningvar_edu0
  outcome = finaldata$hypdrug
  
  outLow <- RDD_educ(Y = outcome)
  
  resL <- outLow$ret
  
  resL <- resL %>%
    rename(Low = res)%>%
    add_rownames

  figL <- c("Low", "1 total",(outLow$estimate)*100, (outLow$CIl)*100, (outLow$CIu)*100)

# Those with high school
  finaldata <- CARRSfull[edu == 1,]

  X = finaldata$runningvar_edu1
  outcome = finaldata$hypdrug
  
  outMiddle <- RDD_educ(Y = outcome)

  resM <- outMiddle$ret
  
  resM <- resM %>%
    rename(Middle = res)

  figM <- c("Middle", "1 total",(outMiddle$estimate)*100, (outMiddle$CIl)*100, (outMiddle$CIu)*100)

  # Those with more than high school 
  finaldata <- CARRSfull[edu == 2,]

  X = finaldata$runningvar_edu2
  outcome = finaldata$hypdrug
  
  outHigh <- RDD_educ(Y = outcome)
  
  resH <- outHigh$ret
  
  resH <- resH %>%
    rename(High = res)

  figH <- c("High", "1 total",(outHigh$estimate)*100, (outHigh$CIl)*100, (outHigh$CIu)*100)

  EduTreat <- cbind(resL,resM,resH)
  figEdu <- rbind(figL, figM, figH)
```

``` {r TreatEduSex, eval = TRUE }
# Up to primary school
  finaldata <- CARRSfull[edu == 0 & female == 0,]

  X = finaldata$runningvar_male_edu0
  outcome = finaldata$hypdrug
  
  outLowM <- RDD_educ(Y = outcome)
  
  resLM <- outLowM$ret
  
  resLM <- resLM %>%
    rename(Low = res)%>%
    add_rownames
  
  figLM <- c("Low", "2 men",(outLowM$estimate)*100, (outLowM$CIl)*100, (outLowM$CIu)*100)  

  
  finaldata <- CARRSfull[edu == 0 & female == 1,]

  X = finaldata$runningvar_female_edu0
  outcome = finaldata$hypdrug
  
  outLowF <- RDD_educ(Y = outcome)
  
  resLF <- outLowF$ret
  
  resLF <- resLF %>%
    rename(Low = res)%>%
    add_rownames

  figLF <- c("Low", "3 women",(outLowF$estimate)*100, (outLowF$CIl)*100, (outLowF$CIu)*100)  


# Those with secondary school
  finaldata <- CARRSfull[edu == 1 & female == 0,]

  X = finaldata$runningvar_male_edu1
  outcome = finaldata$hypdrug
  
  outMiddleM <- RDD_educ(Y = outcome)

  resMM <- outMiddleM$ret
  
  resMM <- resMM %>%
    rename(Middle = res)
 
  figMM <- c("Middle", "2 men",(outMiddleM$estimate)*100, (outMiddleM$CIl)*100, (outMiddleM$CIu)*100)

  
  finaldata <- CARRSfull[edu == 1 & female == 1,]

  X = finaldata$runningvar_female_edu1
  outcome = finaldata$hypdrug
  
  outMiddleF <- RDD_educ(Y = outcome)

  resMF <- outMiddleF$ret
  
  resMF <- resMF %>%
    rename(Middle = res)

  figMF <- c("Middle", "3 women",(outMiddleF$estimate)*100, (outMiddleF$CIl)*100, (outMiddleF$CIu)*100)


  # Those with high school or higher 
  finaldata <- CARRSfull[edu == 2 & female == 0,]

  X = finaldata$runningvar_male_edu2
  outcome = finaldata$hypdrug
  
  outHighM <- RDD_educ(Y = outcome)
  
  resHM <- outHighM$ret
  
  resHM <- resHM %>%
    rename(High = res)

  figHM <- c("High", "2 men",(outHighM$estimate)*100, (outHighM$CIl)*100, (outHighM$CIu)*100)


  
  finaldata <- CARRSfull[edu == 2 & female == 1,]

  X = finaldata$runningvar_female_edu2
  outcome = finaldata$hypdrug
  
  outHighF <- RDD_educ(Y = outcome)
  
  resHF <- outHighF$ret
  
  resHF <- resHF %>%
    rename(High = res)
  
  figHF <- c("High", "3 women",(outHighF$estimate)*100, (outHighF$CIl)*100, (outHighF$CIu)*100)
  
  EduMTreat <- cbind(resLM, resMM, resHM)
  EduFTreat <- cbind(resLF, resMF, resHF)
  
  figSexEdu <- rbind(figLM, figMM, figHM, figLF, figMF, figHF)
  

  # dataset for figures with all results
  TreatmentResultsEdu <- rbind(figEdu,figSexEdu)
  TreatmentResultsEdu <- as.data.frame(TreatmentResultsEdu)
  TreatmentResultsEdu <- rename(TreatmentResultsEdu, edu = V1, bp.sample = V2, estimate = V3, CIl = V4, CIu = V5)
  TreatmentResultsEdu <- transform(TreatmentResultsEdu, estimate = as.numeric(estimate), CIl = as.numeric(CIl), CIu = as.numeric(CIu))
  TreatmentResultsEdu$bp.sample = factor(TreatmentResultsEdu$bp.sample, levels = c("3 women", "2 men", "1 total"))
  TreatmentResultsEdu$edu = factor(TreatmentResultsEdu$edu, levels = c("High", "Middle", "Low"))
```

``` {r TreatAge, eval = TRUE }
# 30-39 years
  finaldata <- CARRSfull[agecat == 0,]

  X = finaldata$runningvar_age0
  outcome = finaldata$hypdrug
  
  outYoung <- RDD_educ(Y = outcome)
  
  resY <- outYoung$ret
  
  resY <- resY %>%
    rename(Young = res)%>%
    add_rownames

  figY <- c("30-39 years", "1 total",(outYoung$estimate)*100, (outYoung$CIl)*100, (outYoung$CIu)*100)

  
# 40+ years
  finaldata <- CARRSfull[agecat == 1,]

  X = finaldata$runningvar_age1
  outcome = finaldata$hypdrug
  
  outOld <- RDD_educ(Y = outcome)

  resO <- outOld$ret
  
  resO <- resO %>%
    rename(Old = res)
  
  figO <- c("40+ years", "1 total",(outOld$estimate)*100, (outOld$CIl)*100, (outOld$CIu)*100)

  AgeTreat <- cbind(resY, resO)
  figAge <- rbind(figY, figO)
```

``` {r TreatAgeSex, eval = TRUE }
# 30-39 years
  finaldata <- CARRSfull[agecat == 0 & female == 0,]

  X = finaldata$runningvar_male_age0
  outcome = finaldata$hypdrug
  
  outYoungM <- RDD_educ(Y = outcome)
  
  resYM <- outYoungM$ret
  
  resYM <- resYM %>%
    rename(Young = res)%>%
    add_rownames
  
  figYM <- c("30-39 years", "2 men",(outYoungM$estimate)*100, (outYoungM$CIl)*100, (outYoungM$CIu)*100)  

    
  finaldata <- CARRSfull[agecat == 0 & female == 1,]

  X = finaldata$runningvar_female_age0
  outcome = finaldata$hypdrug
  
  outYoungF <- RDD_educ(Y = outcome)
  
  resYF <- outYoungF$ret
  
  resYF <- resYF %>%
    rename(Young = res)%>%
    add_rownames
  
  figYF <- c("30-39 years", "3 women",(outYoungF$estimate)*100, (outYoungF$CIl)*100, (outYoungF$CIu)*100)  

# 40+ years
  finaldata <- CARRSfull[agecat == 1 & female == 0,]

  X = finaldata$runningvar_male_age1
  outcome = finaldata$hypdrug
  
  outOldM <- RDD_educ(Y = outcome)

  resOM <- outOldM$ret
  
  resOM <- resOM %>%
    rename(Old = res)

  figOM <- c("40+ years", "2 men",(outOldM$estimate)*100, (outOldM$CIl)*100, (outOldM$CIu)*100)
  
  
  finaldata <- CARRSfull[agecat == 1 & female == 1,]

  X = finaldata$runningvar_female_age1
  outcome = finaldata$hypdrug
  
  outOldF <- RDD_educ(Y = outcome)

  resOF <- outOldF$ret
  
  resOF <- resOF %>%
    rename(Old = res)

  figOF <- c("40+ years", "3 women",(outOldF$estimate)*100, (outOldF$CIl)*100, (outOldF$CIu)*100)
  
  
  AgeMTreat <- cbind(resYM, resOM)
  AgeFTreat <- cbind(resYF, resOF)
  figSexAge <- rbind(figYM, figOM, figYF, figOF)

  # dataset for figures with all results
  TreatmentResultsAge <- rbind(figAge,figSexAge)
  TreatmentResultsAge <- as.data.frame(TreatmentResultsAge)
  TreatmentResultsAge <- rename(TreatmentResultsAge, age = V1, bp.sample = V2, estimate = V3, CIl = V4, CIu = V5)
  TreatmentResultsAge <- transform(TreatmentResultsAge, estimate = as.numeric(estimate), CIl = as.numeric(CIl), CIu = as.numeric(CIu))
  TreatmentResultsAge$bp.sample = factor(TreatmentResultsAge$bp.sample, levels = c("3 women", "2 men", "1 total"))
  TreatmentResultsAge$age = factor(TreatmentResultsAge$age, levels = c("40+ years","30-39 years"))
```

```{r MergingResults, eval = TRUE}
# Education 
  EduDiag$id  <- 1:nrow(EduDiag)
  Edu <- merge(EduDiag,EduTreat, by="rowname")
  Edu <- Edu[order(Edu$id), ]
  Edu <- subset(Edu, select = -c(id))
  
  EduFDiag$id  <- 1:nrow(EduFDiag)
  EduF <- merge(EduFDiag,EduFTreat, by="rowname")
  EduF <- EduF[order(EduF$id), ]
  EduF <- subset(EduF, select = -c(id))
  
  EduMDiag$id  <- 1:nrow(EduMDiag)
  EduM <- merge(EduMDiag,EduMTreat, by="rowname")
  EduM <- EduM[order(EduM$id), ]
  EduM <- subset(EduM, select = -c(id))

# Age  
  AgeDiag$id  <- 1:nrow(AgeDiag)
  Age <- merge(AgeDiag,AgeTreat, by="rowname")
  Age <- Age[order(Age$id), ]
  Age <- subset(Age, select = -c(id))
  
  AgeFDiag$id  <- 1:nrow(AgeFDiag)
  AgeF <- merge(AgeFDiag,AgeFTreat, by="rowname")
  AgeF <- AgeF[order(AgeF$id), ]
  AgeF <- subset(AgeF, select = -c(id))
  
  AgeMDiag$id  <- 1:nrow(AgeMDiag)
  AgeM <- merge(AgeMDiag,AgeMTreat, by="rowname")
  AgeM <- AgeM[order(AgeM$id), ]
  AgeM <- subset(AgeM, select = -c(id))
```

```{r EduTable, eval = TRUE}
#combining diagnosis and treatment results for all three samples
filename <- paste0("edu_combined.tex")
caption=paste0("Effect of home-based screening on hypertension diagnosis and treatment by education categories")

comb_edu_table <- rbind(Edu,EduF, EduM)
rownames(comb_edu_table) <- NULL

combined<-kbl(comb_edu_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccccc", format = "latex",col.names = c("", "Low","Middle","High","Low","Middle","High"),digits = 2)  %>%
    add_header_above(c(" "=1, "Diagnosis" = 3, "Treatment" = 3)) %>%
  kable_styling(latex_options = "HOLD_position", position = "center") %>%
  pack_rows("Panel A: Total", 1,5) %>%
  pack_rows("Estimation results", 1,3) %>%
  pack_rows("Bandwidth details", 4, 5) %>%
  pack_rows(" ", 6,10) %>%
  pack_rows("Panel B: Female", 6,10) %>%
  pack_rows("Estimation results", 6, 8) %>%
  pack_rows("Bandwidth details", 9, 10) %>%
  pack_rows(" ", 11,15) %>%
  pack_rows("Panel C: Male", 11,15) %>%
  pack_rows("Estimation results", 11, 13) %>%
  pack_rows("Bandwidth details", 14, 15) %>%
  footnote(general = c("Note: The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of 0 and confidence intervals resulting from robust bias-corrected standard errors clustered at the individual-level. Low education refers to education up to primary school completed, middle education to secondary school completed, and high education to high school completed or higher.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T) %>%
save_kable(filename, format = "latex")
```

```{r AgeTable, eval = TRUE}
#combining treatment
filename <- paste0("age_combined.tex")
caption=paste0("Effect of home-based screening on hypertension diagnosis and treatment by age group")

comb_age_table <- rbind(Age,AgeF, AgeM)
rownames(comb_age_table) <- NULL

combined<-kbl(comb_age_table,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",col.names = c("","30-39 years","40+ years","30-39 years","40+ years"),digits = 2) %>%
  add_header_above(c(" "=1, "Diagnosis" = 2, "Treatment" = 2)) %>%
  kable_styling(latex_options = "HOLD_position", position = "center") %>%
  pack_rows("Panel A: Total", 1,5) %>%
  pack_rows("Estimation results", 1,3) %>%
  pack_rows("Bandwidth details", 4, 5) %>%
  pack_rows(" ", 6,10) %>%
  pack_rows("Panel B: Female", 6,10) %>%
  pack_rows("Estimation results", 6, 8) %>%
  pack_rows("Bandwidth details", 9, 10) %>%
  pack_rows(" ", 11,15) %>%
  pack_rows("Panel C: Male", 11,15) %>%
  pack_rows("Estimation results", 11, 13) %>%
  pack_rows("Bandwidth details", 14, 15) %>%
  footnote(general = c("Note: The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and confidence intervals resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T) %>%
save_kable(filename, format = "latex")
```


``` {r CoefPlot, eval = TRUE}
# limits of x-axis (originally y axis before flipping)
    ylimedu <- max(max(abs(DiagnosisResultsEdu$CIl)),max(abs(TreatmentResultsEdu$CIl)),
                max(abs(DiagnosisResultsEdu$CIu)), max(abs(TreatmentResultsEdu$CIu)))

# limits of x-axis (originally y axis before flipping)
    ylimage <- max(max(abs(DiagnosisResultsAge$CIl)),max(abs(TreatmentResultsAge$CIl)),
                max(abs(DiagnosisResultsAge$CIu)), max(abs(TreatmentResultsAge$CIu)))

# same limit for education and age plot
    ylim <- ifelse(ylimedu>ylimage,ylimedu,ylimage)

### Education
## Diagnosis
    diaedu <- ggplot(data = DiagnosisResultsEdu, mapping = aes(x = edu, y = estimate, ymin = CIl, ymax = CIu, color = bp.sample)) +
    geom_point(position = position_dodge(0.7), size = 2.5) +
    geom_errorbar(position = position_dodge(0.7), width = 0, size = 0.75) + 
    coord_flip() +
    ylim(-ylim,ylim) +
    geom_hline(yintercept = 0, color="#808080", lwd=1, linetype='dotted') +
    labs(x = "Education categories", y = "Probability (%)\n\n", title = "Diagnosis", family="serif") + 
    theme(text = element_text(size = 24),
          legend.title = element_blank(),
          plot.title.position = "plot",
          plot.title = element_text(size=20,hjust = 0.4,family="serif",face = "bold"),
          axis.title = element_text(size = 20,family="serif"),
          axis.title.y = element_text(margin=margin(r=-100)),
          axis.text.y=element_blank(),
          axis.ticks.y = element_blank(),
          legend.text=element_text(family="serif"),
          legend.position = "none",
          axis.text = element_text(family="serif")) + 
          scale_colour_manual(name = "Population",
                      labels = c("Females", "Males", "Total population"),
                      values = c("#46337EFF", "#B63679FF", "#25AC82FF"),
                      guide = guide_legend(reverse = TRUE)) +   
    # remove grid and grey background
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(),legend.key = element_blank(),
          legend.position="none")


## Treatment
    trtedu <- ggplot(data = TreatmentResultsEdu, mapping = aes(x = edu, y = estimate, ymin = CIl, ymax = CIu, color = bp.sample)) +
    geom_point(position = position_dodge(0.7), size = 2.5) +
    geom_errorbar(position = position_dodge(0.7), width = 0, size = 0.75) + 
    coord_flip() +
    ylim(-ylim,ylim) +
    geom_hline(yintercept = 0, color="#808080", lwd=1, linetype='dotted') + # line at 0  
    labs(x = NULL, y = "Probability (%)\n\n", title = "Treatment", family="serif") + 
    theme(text = element_text(size = 24),
          legend.title = element_blank(),
          plot.title.position = "plot",
          plot.title = element_text(size=20,hjust = 0.5,family="serif",face = "bold"),
          axis.title = element_text(size = 20,family="serif"),
          axis.ticks.y = element_blank(),
          axis.text = element_text(family="serif"),
          axis.text.y = element_text(hjust = 0.5, margin=margin(0,-46,0,0)), # centers labels of y axis and moves panel to align with age plot
          legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
          legend.text=element_text(family="serif"),
          legend.position = "none") +
          scale_colour_manual(name = "Population",
                      labels = c("Females", "Males", "Total population"),
                      values = c("#46337EFF", "#B63679FF", "#25AC82FF"),
                      guide = guide_legend(reverse = TRUE)) +
    # remove grid and grey background
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(),legend.key = element_blank()
          )
educomb <- ggarrange(diaedu, trtedu,
                    ncol = 2, nrow = 1, 
                    widths = c(3.7, 4.5))
educomb

### Age
## Diagnosis
  diaage <- ggplot(data = DiagnosisResultsAge, mapping = aes(x = age, y = estimate, ymin = CIl, ymax = CIu, color = bp.sample)) +
    geom_point(position = position_dodge(0.7), size = 2.5) +
    geom_errorbar(position = position_dodge(0.7), width = 0, size = 0.75) + 
    coord_flip() +
    ylim(-ylim,ylim) +
    geom_hline(yintercept = 0, color="#808080", lwd=1, linetype='dotted') + # line at 0  
    labs(x = "Age groups", y = "Probability (%)",family="serif") + 
    theme(text = element_text(size = 24),
          legend.title = element_blank(),
          plot.title.position = "plot",
          plot.title = element_text(size=20,hjust = 0.5,family="serif"),
          axis.title = element_text(size = 20,family="serif"),
          axis.title.y = element_text(margin=margin(r=-100)),
          axis.text.y=element_blank(),
          axis.ticks.y = element_blank(),
          axis.text = element_text(family="serif"),
          legend.text=element_text(family="serif")) + 
    scale_colour_manual(name = "Population",
                      labels = c("Females", "Males", "Total population"),
                      values = c("#46337EFF", "#B63679FF", "#25AC82FF"),
                      guide = guide_legend(reverse = TRUE)) +   
    # remove grid and grey background
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(),legend.key = element_blank(),
          legend.position="none")

## Treatment  
  trtage <- ggplot(data = TreatmentResultsAge, mapping = aes(x = age, y = estimate, ymin = CIl, ymax = CIu, color = bp.sample)) +
    geom_point(position = position_dodge(0.7), size = 2.5) +
    geom_errorbar(position = position_dodge(0.7), width = 0, size = 0.75) + 
    coord_flip() +
    ylim(-ylim,ylim) +
    geom_hline(yintercept = 0, color="#808080", lwd=1, linetype='dotted') + # line at 0  
    labs(x = NULL, y = "Probability (%)",family="serif") + 
    theme(text = element_text(size = 24),
          legend.title = element_blank(),
          plot.title.position = "plot",
          plot.title = element_text(size=20,hjust = 0.5,family="serif"),
          axis.title = element_text(size = 20,family="serif"),
          axis.ticks.y = element_blank(),
          axis.text.y = element_text(hjust = 0.5, margin=margin(0,-80,0,0)), # centers labels of y axis and moves them a little to the right
          axis.text = element_text(family="serif"),
          legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
          legend.text=element_text(family="serif")) +
    scale_colour_manual(name = "Population",
                      labels = c("Females", "Males", "Total population"),
                      values = c("#46337EFF", "#B63679FF", "#25AC82FF"),
                      guide = guide_legend(reverse = TRUE)) +
    # remove grid and grey background
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(),legend.key = element_blank()
          )

agecomb <- ggarrange(diaage, trtage,
                    ncol = 2, nrow = 1,
                    common.legend = TRUE, legend = "bottom", 
                    widths = c(3.7, 4.5))
agecomb

comb <- ggarrange(educomb, agecomb,
                    nrow = 2,
                    common.legend = TRUE, legend = "bottom",
                    widths = c(3.7, 4.5))

comb
ggsave(paste0("ByEduAge.pdf"), width = 14, height = 10, unit = "in", device = cairo_pdf)  
ggsave(paste0("ByEduAge.jpeg"), width = 14, height = 10, unit = "in", type = cairo)   
```

# RDD function - continuous outcomes

``` {r RDfunctionBP, eval = TRUE }
RDD_educBP <- function(Y) {
  out =  rdrobust(y = Y, X, c = 0, cluster = finaldata$pid, covs = cbind(finaldata$w02, finaldata$w24), 
                  kernel = "triangular", p = 1, rho = 1, all = T)

  ret <- data.frame(
    res = c(formatC(round(out$coef[1,1],2),2,format="f"), formatC(round(out$pv[3,1],2),2,format="f"), 
            paste0( "(", formatC(round(out$ci[3,1],2),2,format="f"), " - ", formatC(round(out$ci[3,2],2),2,format="f"), ")" ), 
            formatC(round(out$bws[1,1],2),2,format="f"), formatC(out$N_h[1] + out$N_h[2]))
  )
  rownames(ret) <- c("Coefficient", "p-value","95\\% CI", "Bandwidth", "Obs. within bandwidth")

  # needed for the figures
  estimate <- out$coef[1,1]
  CIl <- out$ci[3,1]
  CIu <- out$ci[3,2]
  
  return(list(ret=ret, estimate=estimate, CIl=CIl, CIu=CIu))
}
```

# Outcome: systolic blood pressure

``` {r SysEdu, eval = TRUE }
# Those less than secondary school
  finaldata <- CARRSbp[edu == 0,]

  X = finaldata$runningvar_edu0
  outcome = finaldata$sysdiff
  outLow <- RDD_educBP(Y = outcome)
  
  resL <- outLow$ret
  
  resL <- resL %>%
    rename(LowSys = res) %>%
    add_rownames 

  figL <- c("Low", "1 total",(outLow$estimate), (outLow$CIl), (outLow$CIu))
  

# Those with high school
  finaldata <- CARRSbp[edu == 1,]

  X = finaldata$runningvar_edu1
  outcome = finaldata$sysdiff
  
  outMiddle <- RDD_educBP(Y = outcome)

  resM <- outMiddle$ret
  
  resM <- resM %>%
    rename(MiddleSys = res)

  figM <- c("Middle", "1 total",(outMiddle$estimate), (outMiddle$CIl), (outMiddle$CIu))

  # Those with more than high school 
  finaldata <- CARRSbp[edu == 2,]

  X = finaldata$runningvar_edu2
  outcome = finaldata$sysdiff
  
  outHigh <- RDD_educBP(Y = outcome)
  
  resH <- outHigh$ret
  
  resH <- resH %>%
    rename(HighSys = res)

  figH <- c("High", "1 total",(outHigh$estimate), (outHigh$CIl), (outHigh$CIu))

  Edusys <- cbind(resL,resM,resH)
  figEdu <- rbind(figL, figM, figH)
```

``` {r SysEduSex, eval = TRUE }
# Those with primary school or less
  finaldata <- CARRSbp[edu == 0 & female == 0,]

  X = finaldata$runningvar_male_edu0
  outcome = finaldata$sysdiff
  
  outLowM <- RDD_educBP(Y = outcome)
  
  resLM <- outLowM$ret
  
  resLM <- resLM %>%
    rename(LowSys = res)%>%
    add_rownames 

  
  figLM <- c("Low", "2 men",(outLowM$estimate), (outLowM$CIl), (outLowM$CIu))  

  
  finaldata <- CARRSbp[edu == 0 & female == 1,]

  X = finaldata$runningvar_female_edu0
  outcome = finaldata$sysdiff
  
  outLowF <- RDD_educBP(Y = outcome)
  
  resLF <- outLowF$ret
  
  resLF <- resLF %>%
    rename(LowSys = res)%>%
    add_rownames 


  figLF <- c("Low", "3 women",(outLowF$estimate), (outLowF$CIl), (outLowF$CIu))  

# Those with secondary school
  finaldata <- CARRSbp[edu == 1 & female == 0,]

  X = finaldata$runningvar_male_edu1
  outcome = finaldata$sysdiff
  
  outMiddleM <- RDD_educBP(Y = outcome)

  resMM <- outMiddleM$ret
  
  resMM <- resMM %>%
    rename(MiddleSys = res)
 
  figMM <- c("Middle", "2 men",(outMiddleM$estimate), (outMiddleM$CIl), (outMiddleM$CIu))
  
   
  finaldata <- CARRSbp[edu == 1 & female == 1,]

  X = finaldata$runningvar_female_edu1
  outcome = finaldata$sysdiff
  
  outMiddleF <- RDD_educBP(Y = outcome)

  resMF <- outMiddleF$ret
  
  resMF <- resMF %>%
    rename(MiddleSys = res)

  figMF <- c("Middle", "3 women",(outMiddleF$estimate), (outMiddleF$CIl), (outMiddleF$CIu))

  # Those with more than high school 
  finaldata <- CARRSbp[edu == 2 & female == 0,]

  X = finaldata$runningvar_male_edu2
  outcome = finaldata$sysdiff
  
  outHighM <- RDD_educBP(Y = outcome)
  
  resHM <- outHighM$ret
  
  resHM <- resHM %>%
    rename(HighSys = res)

  figHM <- c("High", "2 men",(outHighM$estimate), (outHighM$CIl), (outHighM$CIu))

  
  finaldata <- CARRSbp[edu == 2 & female == 1,]

  X = finaldata$runningvar_female_edu2
  outcome = finaldata$sysdiff
  
  outHighF <- RDD_educBP(Y = outcome)
  
  resHF <- outHighF$ret
  
  resHF <- resHF %>%
    rename(HighSys = res)
  
  figHF <- c("High", "3 women",(outHighF$estimate), (outHighF$CIl), (outHighF$CIu))
  
  EduMsys <- cbind(resLM, resMM, resHM)
  EduFsys <- cbind(resLF, resMF, resHF)

  figSexEdu <- rbind(figLM, figMM, figHM, figLF, figMF, figHF)
  

  # dataset for figures with all results
  SystolicResultsEdu <- rbind(figEdu,figSexEdu)
  SystolicResultsEdu <- as.data.frame(SystolicResultsEdu)
  SystolicResultsEdu <- rename(SystolicResultsEdu, edu = V1, bp.sample = V2, estimate = V3, CIl = V4, CIu = V5)
  SystolicResultsEdu <- transform(SystolicResultsEdu, estimate = as.numeric(estimate), CIl = as.numeric(CIl), CIu = as.numeric(CIu))
  SystolicResultsEdu$bp.sample = factor(SystolicResultsEdu$bp.sample, levels = c("3 women", "2 men", "1 total"))
  SystolicResultsEdu$edu = factor(SystolicResultsEdu$edu, levels = c("High", "Middle", "Low"))
```

``` {r SysAge, eval = TRUE }
# 30-39 years
  finaldata <- CARRSbp[agecat == 0,]

  X = finaldata$runningvar_age0
  outcome = finaldata$sysdiff
  outYoung <- RDD_educBP(Y = outcome)
  
  resY <- outYoung$ret
  
  resY <- resY %>%
    rename(Young = res)%>%
    add_rownames 

  figY <- c("30-39 years", "1 total",(outYoung$estimate), (outYoung$CIl), (outYoung$CIu))

# 40+ years
  finaldata <- CARRSbp[agecat == 1,]

  X = finaldata$runningvar_age1
  outcome = finaldata$sysdiff
  
  outOld <- RDD_educBP(Y = outcome)

  resO <- outOld$ret
  
  resO <- resO %>%
    rename(Old = res)
  
  figO <- c("40+ years", "1 total",(outOld$estimate), (outOld$CIl), (outOld$CIu))

  Agesys <- cbind(resY, resO)
  figAge <- rbind(figY, figO)
```

``` {r SysAgeSex, eval = TRUE }
# 30-39 years
finaldata <- CARRSbp[agecat == 0 & female == 0,]

  X = finaldata$runningvar_male_age0
  outcome = finaldata$sysdiff
  
  outYoungM <- RDD_educBP(Y = outcome)
  
  resYM <- outYoungM$ret
  
  resYM <- resYM %>%
    rename(Young = res)%>%
    add_rownames 
  
  figYM <- c("30-39 years", "2 men",(outYoungM$estimate), (outYoungM$CIl), (outYoungM$CIu))  

finaldata <- CARRSbp[agecat == 0 & female == 1,]

  X = finaldata$runningvar_female_age0
  outcome = finaldata$sysdiff
  
  outYoungF <- RDD_educBP(Y = outcome)
  
  resYF <- outYoungF$ret
  
  resYF <- resYF %>%
    rename(Young = res)%>%
    add_rownames 

  figYF <- c("30-39 years", "3 women",(outYoungF$estimate), (outYoungF$CIl), (outYoungF$CIu))  

# 40+ years
  finaldata <- CARRSbp[agecat == 1 & female == 0,]

  X = finaldata$runningvar_male_age1
  outcome = finaldata$sysdiff
  
  outOldM <- RDD_educBP(Y = outcome)

  resOM <- outOldM$ret
  
  resOM <- resOM %>%
    rename(Old = res)

  figOM <- c("40+ years", "2 men",(outOldM$estimate), (outOldM$CIl), (outOldM$CIu))
  
  finaldata <- CARRSbp[agecat == 1 & female == 1,]

  X = finaldata$runningvar_female_age1
  outcome = finaldata$sysdiff
  
  outOldF <- RDD_educBP(Y = outcome)

  resOF <- outOldF$ret
  
  resOF <- resOF %>%
    rename(Old = res)

  figOF <- c("40+ years", "3 women",(outOldF$estimate), (outOldF$CIl), (outOldF$CIu))
  
  AgeMsys <- cbind(resYM, resOM)
  AgeFsys <- cbind(resYF, resOF)
  
  figSexAge <- rbind(figYM, figOM, figYF, figOF)

  # dataset for figures with all results
  SystolicResultsAge <- rbind(figAge,figSexAge)
  SystolicResultsAge <- as.data.frame(SystolicResultsAge)
  SystolicResultsAge <- rename(SystolicResultsAge, age = V1, bp.sample = V2, estimate = V3, CIl = V4, CIu = V5)
  SystolicResultsAge <- transform(SystolicResultsAge, estimate = as.numeric(estimate), CIl = as.numeric(CIl), CIu = as.numeric(CIu))
  SystolicResultsAge$bp.sample = factor(SystolicResultsAge$bp.sample, levels = c("3 women", "2 men", "1 total"))
  SystolicResultsAge$age = factor(SystolicResultsAge$age, levels = c("40+ years","30-39 years"))
```

# Outcome: diastolic blood pressure

``` {r DiasEdu, eval = TRUE }
# Those less than secondary school
  finaldata <- CARRSbp[edu == 0,]

  X = finaldata$runningvar_edu0
  outcome = finaldata$diasdiff
  
  outLow <- RDD_educBP(Y = outcome)
  
  resL <- outLow$ret
  
  resL <- resL %>%
    rename(Low = res)%>%
    add_rownames 


  figL <- c("Low", "1 total",(outLow$estimate), (outLow$CIl), (outLow$CIu))

# Those with high school
  finaldata <- CARRSbp[edu == 1,]

  X = finaldata$runningvar_edu1
  outcome = finaldata$diasdiff
  
  outMiddle <- RDD_educBP(Y = outcome)

  resM <- outMiddle$ret
  
  resM <- resM %>%
    rename(Middle = res)

  figM <- c("Middle", "1 total",(outMiddle$estimate), (outMiddle$CIl), (outMiddle$CIu))

  # Those with more than high school 
  finaldata <- CARRSbp[edu == 2,]

  X = finaldata$runningvar_edu2
  outcome = finaldata$diasdiff
  
  outHigh <- RDD_educBP(Y = outcome)
  
  resH <- outHigh$ret
  
  resH <- resH %>%
    rename(High = res)

  figH <- c("High", "1 total",(outHigh$estimate), (outHigh$CIl), (outHigh$CIu))

  Edudias <- cbind(resL,resM,resH)
  figEdu <- rbind(figL, figM, figH)
```

``` {r DiasEduSex, eval = TRUE }
# Up to primary school
  finaldata <- CARRSbp[edu == 0 & female == 0,]

  X = finaldata$runningvar_male_edu0
  outcome = finaldata$diasdiff
  
  outLowM <- RDD_educBP(Y = outcome)
  
  resLM <- outLowM$ret
  
  resLM <- resLM %>%
    rename(Low = res)%>%
    add_rownames 

  
  figLM <- c("Low", "2 men",(outLowM$estimate), (outLowM$CIl), (outLowM$CIu))  

    
  finaldata <- CARRSbp[edu == 0 & female == 1,]

  X = finaldata$runningvar_female_edu0
  outcome = finaldata$diasdiff
  
  outLowF <- RDD_educBP(Y = outcome)
  
  resLF <- outLowF$ret
  
  resLF <- resLF %>%
    rename(Low = res)%>%
    add_rownames 


  figLF <- c("Low", "3 women",(outLowF$estimate), (outLowF$CIl), (outLowF$CIu))  


# Those with secondary school
  finaldata <- CARRSbp[edu == 1 & female == 0,]

  X = finaldata$runningvar_male_edu1
  outcome = finaldata$diasdiff
  
  outMiddleM <- RDD_educBP(Y = outcome)

  resMM <- outMiddleM$ret
  
  resMM <- resMM %>%
    rename(Middle = res)
 
  figMM <- c("Middle", "2 men",(outMiddleM$estimate), (outMiddleM$CIl), (outMiddleM$CIu))

  
  finaldata <- CARRSbp[edu == 1 & female == 1,]

  X = finaldata$runningvar_female_edu1
  outcome = finaldata$diasdiff
  
  outMiddleF <- RDD_educBP(Y = outcome)

  resMF <- outMiddleF$ret
  
  resMF <- resMF %>%
    rename(Middle = res)

  figMF <- c("Middle", "3 women",(outMiddleF$estimate), (outMiddleF$CIl), (outMiddleF$CIu))


  # Those with high school or higher 
  finaldata <- CARRSbp[edu == 2 & female == 0,]

  X = finaldata$runningvar_male_edu2
  outcome = finaldata$diasdiff
  
  outHighM <- RDD_educBP(Y = outcome)
  
  resHM <- outHighM$ret
  
  resHM <- resHM %>%
    rename(High = res)

  figHM <- c("High", "2 men",(outHighM$estimate), (outHighM$CIl), (outHighM$CIu))


  
  finaldata <- CARRSbp[edu == 2 & female == 1,]

  X = finaldata$runningvar_female_edu2
  outcome = finaldata$diasdiff
  
  outHighF <- RDD_educBP(Y = outcome)
  
  resHF <- outHighF$ret
  
  resHF <- resHF %>%
    rename(High = res)
  
  figHF <- c("High", "3 women",(outHighF$estimate), (outHighF$CIl), (outHighF$CIu))
  
  EduMdias <- cbind(resLM, resMM, resHM)
  EduFdias <- cbind(resLF, resMF, resHF)
  
  figSexEdu <- rbind(figLM, figMM, figHM, figLF, figMF, figHF)

  # dataset for figures with all results
  DiastolicResultsEdu <- rbind(figEdu,figSexEdu)
  DiastolicResultsEdu <- as.data.frame(DiastolicResultsEdu)
  DiastolicResultsEdu <- rename(DiastolicResultsEdu, edu = V1, bp.sample = V2, estimate = V3, CIl = V4, CIu = V5)
  DiastolicResultsEdu <- transform(DiastolicResultsEdu, estimate = as.numeric(estimate), CIl = as.numeric(CIl), CIu = as.numeric(CIu))
  DiastolicResultsEdu$bp.sample = factor(DiastolicResultsEdu$bp.sample, levels = c("3 women", "2 men", "1 total"))
  DiastolicResultsEdu$edu = factor(DiastolicResultsEdu$edu, levels = c("High", "Middle", "Low"))
```

``` {r DiasAge, eval = TRUE }
# 30-39 years
  finaldata <- CARRSbp[agecat == 0,]

  X = finaldata$runningvar_age0
  outcome = finaldata$diasdiff
  
  outYoung <- RDD_educBP(Y = outcome)
  
  resY <- outYoung$ret
  
  resY <- resY %>%
    rename(Young = res)%>%
    add_rownames 

  figY <- c("30-39 years", "1 total",(outYoung$estimate), (outYoung$CIl), (outYoung$CIu))

  
# 40+ years
  finaldata <- CARRSbp[agecat == 1,]

  X = finaldata$runningvar_age1
  outcome = finaldata$diasdiff
  
  outOld <- RDD_educBP(Y = outcome)

  resO <- outOld$ret
  
  resO <- resO %>%
    rename(Old = res)
  
  figO <- c("40+ years", "1 total",(outOld$estimate), (outOld$CIl), (outOld$CIu))

  Agedias <- cbind(resY, resO)
  figAge <- rbind(figY, figO)
```

``` {r DiasAgeSex, eval = TRUE }
# 30-39 years
  finaldata <- CARRSbp[agecat == 0 & female == 0,]

  X = finaldata$runningvar_male_age0
  outcome = finaldata$diasdiff
  
  outYoungM <- RDD_educBP(Y = outcome)
  
  resYM <- outYoungM$ret
  
  resYM <- resYM %>%
    rename(Young = res)%>%
    add_rownames 
  
  figYM <- c("30-39 years", "2 men",(outYoungM$estimate), (outYoungM$CIl), (outYoungM$CIu))  

    
  finaldata <- CARRSbp[agecat == 0 & female == 1,]

  X = finaldata$runningvar_female_age0
  outcome = finaldata$diasdiff
  
  outYoungF <- RDD_educBP(Y = outcome)
  
  resYF <- outYoungF$ret
  
  resYF <- resYF %>%
    rename(Young = res)%>%
    add_rownames 

  figYF <- c("30-39 years", "3 women",(outYoungF$estimate), (outYoungF$CIl), (outYoungF$CIu))  


# 40+ years
  finaldata <- CARRSbp[agecat == 1 & female == 0,]

  X = finaldata$runningvar_male_age1
  outcome = finaldata$diasdiff
  
  outOldM <- RDD_educBP(Y = outcome)

  resOM <- outOldM$ret
  
  resOM <- resOM %>%
    rename(Old = res)

  figOM <- c("40+ years", "2 men",(outOldM$estimate), (outOldM$CIl), (outOldM$CIu))
  
  
  finaldata <- CARRSbp[agecat == 1 & female == 1,]

  X = finaldata$runningvar_female_age1
  outcome = finaldata$diasdiff
  
  outOldF <- RDD_educBP(Y = outcome)

  resOF <- outOldF$ret
  
  resOF <- resOF %>%
    rename(Old = res)

  figOF <- c("40+ years", "3 women",(outOldF$estimate), (outOldF$CIl), (outOldF$CIu))
  
  
  AgeMdias <- cbind(resYM, resOM)
  AgeFdias <- cbind(resYF, resOF)
  
  figSexAge <- rbind(figYM, figOM, figYF, figOF)

  # dataset for figures with all results
  DiastolicResultsAge <- rbind(figAge,figSexAge)
  DiastolicResultsAge <- as.data.frame(DiastolicResultsAge)
  DiastolicResultsAge <- rename(DiastolicResultsAge, age = V1, bp.sample = V2, estimate = V3, CIl = V4, CIu = V5)
  DiastolicResultsAge <- transform(DiastolicResultsAge, estimate = as.numeric(estimate), CIl = as.numeric(CIl), CIu = as.numeric(CIu))
  DiastolicResultsAge$bp.sample = factor(DiastolicResultsAge$bp.sample, levels = c("3 women", "2 men", "1 total"))
  DiastolicResultsAge$age = factor(DiastolicResultsAge$age, levels = c("40+ years","30-39 years"))
```


``` {r MergingResultsBP, eval = TRUE}
# Education 
  Edusys$id  <- 1:nrow(Edusys)
  Edu_bp <- merge(Edusys,Edudias, by="rowname")
  Edu_bp <- Edu_bp[order(Edu_bp$id), ]
  Edu_bp <- subset(Edu_bp, select = -c(id))
  
  EduFsys$id  <- 1:nrow(EduFsys)
  EduF_bp <- merge(EduFsys,EduFdias, by="rowname")
  EduF_bp <- EduF_bp[order(EduF_bp$id), ]
  EduF_bp <- subset(EduF_bp, select = -c(id))
  
  EduMsys$id  <- 1:nrow(EduMsys)
  EduM_bp <- merge(EduMsys,EduMdias, by="rowname")
  EduM_bp <- EduM_bp[order(EduM_bp$id), ]
  EduM_bp <- subset(EduM_bp, select = -c(id))

# Age  
  Agesys$id  <- 1:nrow(Agesys)
  Age_bp <- merge(Agesys,Agedias, by="rowname")
  Age_bp <- Age_bp[order(Age_bp$id), ]
  Age_bp <- subset(Age_bp, select = -c(id))
  
  AgeFsys$id  <- 1:nrow(AgeFsys)
  AgeF_bp <- merge(AgeFsys,AgeFdias, by="rowname")
  AgeF_bp <- AgeF_bp[order(AgeF_bp$id), ]
  AgeF_bp <- subset(AgeF_bp, select = -c(id))
  
  AgeMsys$id  <- 1:nrow(AgeMsys)
  AgeM_bp <- merge(AgeMsys,AgeMdias, by="rowname")
  AgeM_bp <- AgeM_bp[order(AgeM_bp$id), ]
  AgeM_bp <- subset(AgeM_bp, select = -c(id))
```


``` {r EduTableBP, eval = TRUE }
#combining diagnosis and treatment results for all three samples
filename <- paste0("edu_combined_bpdiff.tex")
caption=paste0("Effect of home-based screening on systolic and diastolic blood pressure by education categories")

comb_edu_table_bp <- rbind(Edu_bp,EduF_bp, EduM_bp)
rownames(comb_edu_table_bp) <- NULL

combined<-kbl(comb_edu_table_bp,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccccc", format = "latex",col.names = c("", "Low","Middle","High","Low","Middle","High"),digits = 2)  %>%
    add_header_above(c(" "=1, "Systolic" = 3, "Diastolic" = 3)) %>%
  kable_styling(latex_options = "HOLD_position", position = "center") %>%
  pack_rows("Panel A: Total", 1,5) %>%
  pack_rows("Estimation results", 1,3) %>%
  pack_rows("Bandwidth details", 4, 5) %>%
  pack_rows(" ", 6,10) %>%
  pack_rows("Panel B: Female", 6,10) %>%
  pack_rows("Estimation results", 6, 8) %>%
  pack_rows("Bandwidth details", 9, 10) %>%
  pack_rows(" ", 11,15) %>%
  pack_rows("Panel C: Male", 11,15) %>%
  pack_rows("Estimation results", 11, 13) %>%
  pack_rows("Bandwidth details", 14, 15) %>%
  footnote(general = c("Note: The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of 0 and confidence intervals resulting from robust bias-corrected standard errors clustered at the individual-level. Low education refers to education up to primary school completed, middle education to secondary school completed, and high education to high school completed or higher.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T) %>%
save_kable(filename, format = "latex")
```


``` {r AgeTableBP, eval = TRUE}
filename <- paste0("age_combined_bpdiff.tex")
caption=paste0("Effect of home-based screening on systolic and diastolic blood pressure by age group")

comb_age_table_bp <- rbind(Age_bp,AgeF_bp, AgeM_bp)
rownames(comb_age_table_bp) <- NULL


combined<-kbl(comb_age_table_bp,caption=caption, booktabs = T, linesep = '', escape = FALSE,align = "lccccc", format = "latex",col.names = c("","30-39 years","40+ years","30-39 years","40+ years"),digits = 2) %>%
  add_header_above(c(" "=1, "Systolic" = 2, "Diastolic" = 2)) %>%
  kable_styling(latex_options = "HOLD_position", position = "center") %>%
  pack_rows("Panel A: Total", 1,5) %>%
  pack_rows("Estimation results", 1,3) %>%
  pack_rows("Bandwidth details", 4, 5) %>%
  pack_rows(" ", 6,10) %>%
  pack_rows("Panel B: Female", 6,10) %>%
  pack_rows("Estimation results", 6, 8) %>%
  pack_rows("Bandwidth details", 9, 10) %>%
  pack_rows(" ", 11,15) %>%
  pack_rows("Panel C: Male", 11,15) %>%
  pack_rows("Estimation results", 11, 13) %>%
  pack_rows("Bandwidth details", 14, 15) %>%
  footnote(general = c("Note: The model estimation included a local linear trend, triangular Kernel weights, and a data-driven bandwidth selection (mean squared error optimal bandwidth). The table displays the local average treatment effect in percentage points for individuals with a standardized blood pressure of zero and confidence intervals resulting from robust bias-corrected standard errors clustered at the individual-level.","Abbreviation: CI = confidence interval"),general_title="",threeparttable = T, fixed_small_size = T) %>%
save_kable(filename, format = "latex")

```


``` {r CoefPlotBP, eval = TRUE}
# limits of x-axis (originally y axis before flipping)
    ylimedu <- max(max(abs(SystolicResultsEdu$CIl)),max(abs(DiastolicResultsEdu$CIl)),
                max(abs(SystolicResultsEdu$CIu)), max(abs(DiastolicResultsEdu$CIu)))

# limits of x-axis (originally y axis before flipping)
    ylimage <- max(max(abs(SystolicResultsAge$CIl)),max(abs(DiastolicResultsAge$CIl)),
                max(abs(SystolicResultsAge$CIu)), max(abs(DiastolicResultsAge$CIu)))

# same limit for education and age plot
    ylim <- ifelse(ylimedu>ylimage,ylimedu,ylimage)

### Education
## Systolic BP
    sysedu <- ggplot(data = SystolicResultsEdu, mapping = aes(x = edu, y = estimate, ymin = CIl, ymax = CIu, color = bp.sample)) +
    geom_point(position = position_dodge(0.7), size = 2.5) +
    geom_errorbar(position = position_dodge(0.7), width = 0, size = 0.75) + 
    coord_flip() +
    ylim(-ylim,ylim) +
    geom_hline(yintercept = 0, color="#808080", lwd=1, linetype='dotted') + # line at 0  
    labs(x = "Education categories", y = "Change (mmHg)\n\n", title = "Systolic", family="serif") + 
    theme(text = element_text(size = 24),
          legend.title = element_blank(),
          plot.title.position = "plot",
          plot.title = element_text(size=20,hjust = 0.53,family="serif",face = "bold"),
          plot.margin = margin(r=15), # x-axis labels were cut off
          axis.title = element_text(size = 20,family="serif"),
          axis.text.y=element_blank(),
          axis.ticks.y = element_blank(),
          legend.text=element_text(family="serif"),
          legend.position = "none",
          axis.text = element_text(family="serif")) + 
          scale_colour_manual(name = "Population",
                      labels = c("Females", "Males", "Total population"),
                      values = c("#46337EFF", "#B63679FF", "#25AC82FF"),
                      guide = guide_legend(reverse = TRUE)) +   
    # remove grid and grey background
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(),legend.key = element_blank(),
          legend.position="none")

## Diastolic BP
    diasedu <- ggplot(data = DiastolicResultsEdu, mapping = aes(x = edu, y = estimate, ymin = CIl, ymax = CIu, color = bp.sample)) +
    geom_point(position = position_dodge(0.7), size = 2.5) +
    geom_errorbar(position = position_dodge(0.7), width = 0, size = 0.75) + 
    coord_flip() +
    ylim(-ylim,ylim) +
    geom_hline(yintercept = 0, color="#808080", lwd=1, linetype='dotted') + # line at 0  
    labs(x = NULL, y = "Change (mmHg)\n\n", title = "Diastolic", family="serif") + 
    theme(text = element_text(size = 24),
          legend.title = element_blank(),
          plot.title.position = "plot",
          plot.title = element_text(size=20,hjust = 0.6,family="serif",face = "bold"),
          plot.margin = margin(r=15), # x-axis labels were cut off
          axis.title = element_text(size = 20,family="serif"),
          axis.ticks.y = element_blank(),
          axis.text = element_text(family="serif"),
          axis.text.y = element_text(hjust = 0.5, margin=margin(0,38,0,0)), # centers labels of y axis and moves panel to align with age plot
          legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
          legend.text=element_text(family="serif"),
          legend.position = "none") +
          scale_colour_manual(name = "Population",
                      labels = c("Females", "Males", "Total population"),
                      values = c("#46337EFF", "#B63679FF", "#25AC82FF"),
                      guide = guide_legend(reverse = TRUE)) +
    # remove grid and grey background
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(),legend.key = element_blank()
          )

educombBP <- ggarrange(sysedu, diasedu,
                    ncol = 2, nrow = 1, 
                    widths = c(3.7, 4.5))
educombBP


### Age
## systolic
  sysage <- ggplot(data = SystolicResultsAge, mapping = aes(x = age, y = estimate, ymin = CIl, ymax = CIu, color = bp.sample)) +
    geom_point(position = position_dodge(0.7), size = 2.5) +
    geom_errorbar(position = position_dodge(0.7), width = 0, size = 0.75) + 
    coord_flip() +
    ylim(-ylim,ylim) +
    geom_hline(yintercept = 0, color="#808080", lwd=1, linetype='dotted') + # line at 0  
    labs(x = "Age groups", y = "Change (mmHg)",family="serif") + 
    theme(text = element_text(size = 24),
          legend.title = element_blank(),
          plot.title.position = "plot",
          plot.title = element_text(size=20,hjust = 0.5,family="serif"),
          plot.margin = margin(r=15), # x-axis labels were cut off
          axis.title = element_text(size = 20,family="serif"),
          axis.text.y=element_blank(),
          axis.ticks.y = element_blank(),
          axis.text = element_text(family="serif"),
          legend.text=element_text(family="serif")) + 
    scale_colour_manual(name = "Population",
                      labels = c("Females", "Males", "Total population"),
                      values = c("#46337EFF", "#B63679FF", "#25AC82FF"),
                      guide = guide_legend(reverse = TRUE)) +   
    # remove grid and grey background
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(),legend.key = element_blank(),
          legend.position="none")


## diastolic
  diasage <- ggplot(data = DiastolicResultsAge, mapping = aes(x = age, y = estimate, ymin = CIl, ymax = CIu, color = bp.sample)) +
    geom_point(position = position_dodge(0.7), size = 2.5) +
    geom_errorbar(position = position_dodge(0.7), width = 0, size = 0.75) + 
    coord_flip() +
    ylim(-ylim,ylim) +
    geom_hline(yintercept = 0, color="#808080", lwd=1, linetype='dotted') + # line at 0  
    labs(x = NULL, y = "Change (mmHg)",family="serif") + 
    theme(text = element_text(size = 24),
          legend.title = element_blank(),
          plot.title.position = "plot",
          plot.title = element_text(size=20,hjust = 0.5,family="serif"),
          plot.margin = margin(r=15), # x-axis labels were cut off
          axis.title = element_text(size = 20,family="serif"),
          axis.ticks.y = element_blank(),
          axis.text = element_text(family="serif"),
          legend.background = element_rect(colour = 'black', fill = 'white', linetype='solid'),
          legend.text=element_text(family="serif")) +
    scale_colour_manual(name = "Population",
                      labels = c("Females", "Males", "Total population"),
                      values = c("#46337EFF", "#B63679FF", "#25AC82FF"),
                      guide = guide_legend(reverse = TRUE)) +
    # remove grid and grey background
    theme(panel.grid.major = element_blank(), panel.grid.minor = element_blank(),
          panel.background = element_blank(),legend.key = element_blank()
          )

agecombBP <- ggarrange(sysage, diasage,
                    ncol = 2, nrow = 1,
                    common.legend = TRUE, legend = "bottom", 
                    widths = c(3.7, 4.5))
agecombBP

combBP <- ggarrange(educombBP, agecombBP,
                    nrow = 2,
                    common.legend = TRUE, legend = "bottom",
                    widths = c(3.7, 4.5))

ggsave(paste0("ByEduAge_bpdiff.pdf"), width = 14, height = 10, unit = "in", device = cairo_pdf)  
ggsave(paste0("ByEduAge_bpdiff.jpeg"), width = 14, height = 10, unit = "in", type = cairo)   
```
